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Abstract — Distributed Space Missions (DSMs) are gaining 
momentum in their application to earth science missions owing 
to their unique ability to increase observation sampling in 
spatial, spectral and temporal dimensions simultaneously. 
DSM architectures have a large number of design variables 
and since they are expected to increase mission flexibility, 
scalability, evolvability and robustness, their design is a 
complex problem with many variables and objectives affecting 
performance. There are very few open-access tools available to 
explore the tradespace of variables which allow performance 
assessment and are easy to plug into science goals, and 
therefore select the most optimal design. This paper presents a 
software tool developed on the MATLAB engine interfacing 
with STK, for DSM orbit design and selection. It is capable of 
generating thousands of homogeneous constellation or 
formation flight architectures based on pre-defined design 
variable ranges and sizing those architectures in terms of pre- 
defined performance metrics. The metrics can be input into 
observing system simulation experiments, as available from the 
science teams, allowing dynamic coupling of science and 
engineering designs. Design variables include but are not 
restricted to constellation type, formation flight type, FOV of 
instrument, altitude and inclination of chief orbits, differential 
orbital elements, leader satellites, latitudes or regions of 
interest, planes and satellite numbers. Intermediate 
performance metrics include angular coverage, number of 
accesses, revisit coverage, access deterioration over time at 
every point of the Earth’s grid. The orbit design process can be 
streamlined and variables more bounded along the way, owing 
to the availability of low fidelity and low complexity models 
such as corrected HCW equations up to high precision STK 
models with J2 and drag. The tool can thus help any scientist 
or program manager select pre-Phase A, Pareto optimal DSM 
designs for a variety of science goals without having to delve 
into the details of the engineering design process. 
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1. Introduction 

Distributed Space Missions (DSMs) are gaining momentum 
in their application to Earth Observation (EO) missions 
owing to their unique ability to increase observation 
sampling in spatial, spectral, angular and temporal 
dimensions simultaneously. Spatial resolution of an image 
can be increased by using multiple satellites in formation 
flight to synthesise a long baseline aperture as shown for 
optical interferometry[l] and synthetic aperture radars. 
Constellations of evenly spaced satellites on repeat track 
orbits ensure temporal sampling within a few hours as well 
as continuous coverage maintenance. Spectral sampling can 
be improved by fractionating the payload (fractionated 
spacecraft) such that each physical entity images a different 
part of the spectrum and has customized optics to do so. 
Angular sampling or the ability to look at the same point on 
the ground at different angles (for reflectance studies or 
navigation) improves by flying many satellites in 
formation[2]. Since DSMs allow sampling improvement in 
any dimension by increasing number of satellites instead of 
individual sizes, radiometric resolution can be improved 
without compromising on other sampling requirements. 

DSMs can be considered homogeneous or heterogeneous 
combinations of monoliths. They include homogenous 
constellations like the Global Positioning System or 
heterogeneous ad-hoc flyers like the A-Train, autonomous 
formation flying clusters such as PRISMA, fractionated 
spacecraft such as the recently canceled System F6 
Program[3] and cellularized systems such as the DARPA 
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Phoenix Program[4]. Formation flight, as required in 
clusters, fractionation or cellularization, entails active 
control of the individual spacecrafts in order to maintain 
relative distances, orientations and geometry [5]. 
Fractionated spacecraft have the different spacecraft 
subsystems distributed over the physical entities and they 
exchange data, power and telemetry among each other. 
Cellularized systems are formed by assembling on-orbit 
resources called satlets to make aggregated, distributed 
systems. 

Since DSM architectures are defined by monolithic 
architecture variables and variables associated with the 
distributed framework, it leads to a large number of design 
variables. The number increases further in heterogeneous 
cases. DSMs are also expected to increase mission 
flexibility, scalability, evolvability and robustness as well as 
to minimize costs and risks associated with launch and 
operations. Thus, DSM design is a complex problem with 
many design variables, multiple objectives for performance, 
cost and emergent, often unexpected, behaviors. 

2. Decisions in DSM Design 

Designing space systems is not only technically challenging 
but also involves making hundreds of decisions early in the 
design cycle for allocating limited resources across the 
system and optimizing performance and cost. EO 
performance can be simplistically represented by spatial 
resolution, spatial range (swath, coverage), spectral 
resolution (wavelength bandwidth), spectral range (spectrum 
covered), angular resolution (number of view and solar 
illumination angles for the same image), angular range 
(spread of those angles), temporal range (mission lifetime), 
temporal resolution (repeat or revisit time), radiometric 
range (number of bits) and radiometric resolution (bits, 
signal to noise ratio). 

Distributed systems have all the trades associated with 
monolithic systems and more associated with the network. 
Extra design variables include but are not restricted to the 
number of satellites and their individual masses, their orbits 
and inter- satellite spacing, existence and nature of inter- 
satellite communication and downlink schedules. These 
variables directly impact performance and cost. 
Performance variables, as defined, can be mutually 
conflicting across the spatial, spectral, temporal, angular and 
radiometric dimensions and within each dimension. For 
example, most earth observation satellites are placed in 
repeat ground track orbits so that the same point on the 
Earth is revisited regularly and frequently. Obviously, more 
frequent revisits imply that the rotation of the Earth and the 


orbit has to be adjusted in such a way that the satellite 
comes back to the same spot frequently, and as a result has 
less time to visit similar spots on other longitudes. 
Therefore, global spatial coverage or spatial range and 
temporal resolution are conflicting metrics. Both can be 
improved by increasing the swath of the satellite or the size 
of any instantaneous ground image. However, for a given 
number of pixels in an image, increasing its size or swath 
with increases the size of the pixels and coarsens resolution. 
Therefore, spatial range and temporal resolution are both 
conflicting metrics with respect to spatial resolution. Design 
variables need to be permuted to consider architectures that 
trade these metrics for an optimal design. 

There are very few open-access tools available to explore 
the tradespace of variables, minimize cost and maximize 
performance for pre-defmed science goals, and therefore 
select the most optimal design. This paper concentrates on a 
tool for computing angular performance metrics and their 
trade-offs with spatial and temporal performance. The 
spatial metrics obviously affect the spectral and radiometric 
performance, depending on the kind of spectrometer or 
radiometer used, which has been discussed in previous 
literature [6]. 

3. Distributed spacecraft for multi- 
angular OBSERVATION 

Angular sampling implies taking images of the same ground 
spot at multiple 3D angles of solar incidence and reflection 
simultaneously. A near- simultaneous measurement 
requirement deems monoliths insufficient for accurate and 
dense angular sampling[7]. Monolithic spacecraft have 
traditionally approximated the angular samples by 
combining measurements taken over time with forward-aft 
(e.g. TERRA’S MISR[8]) or cross-track swath (e.g. 
TERRA’S MODIS[9]) sensors. However, a single satellite 
can make measurements only along a restrictive plane with 
respect to the solar phase and most earth observation 
satellites are even more restricted since they are on sun- 
synchronous orbits. Further, the angular measurements are 
separated in time by many minutes along-track or weeks 
cross-track. In areas of fast changing surface/cloud 
conditions especially during the snow melt season/tropical 
storms, a few days can make a big difference in reflectance. 

Near-simultaneous angular sampling can be improved by 
using a cluster or constellation of nanosatellites on a 
repeating-ground-track orbit[7]. The cluster can make multi- 
spectral measurements of a ground spot at multiple 3D 
angles at the same time as they pass overhead either using 
narrow field of view (NFOV) instruments in controlled 
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formation flight (Figure 1-a) or wide field of view (WFOV) 
instruments with overlapping ground spots providing 
integrated images at various angles (Figure 1-b). Parallel 
studies have demonstrated the technical feasibility of 
subsystems [10], suitability of payload development [6] to 
support such a mission, availability of science models to 
quantify the performance of such DSMs[2],[12],[13] as well 
as open-source flight software to continually update satellite 
capability for staged, scalable deployment [14]-[17]. This 
paper focuses on generation of feasible formation flight and 
constellation architectures for all types of multi-angular 
measurements on the earth (Figure 1). 


Figure 1: (a) A DSM making multi-angular, multi- 
spectral measurements by virtue of pointing its NFOVs 
at the same ground spot, as it orbits the Earth as a single 
system (adapted from Leonardo BRDF[18]). (b) A DSM 
making multi-angular, multi-spectral measurements by 
virtue of their overlapping WFOVs at different angles. 

Metrics for Angular Sampling 

The widely accepted metric to quantify the angular 
dependence of remotely sensed signal is called BRDF or 
Bidirectional Reflectance-distribution function. BRDF of an 
optically thick body is a property of the surface material and 
its roughness. It is the ratio of reflected radiance to incident 
irradiance that depends on 3D geometry of incident and 
reflected elementary beams [19]. It depends on four major 
angles - the solar zenith (SZA) and azimuth angle (SAZ) 
and the view zenith (VZA) and azimuth angle (VAZ). The 
azimuth angles are simplified to one angle called the relative 
azimuth angle (RAZ). BRDF is used for the derivation of 
surface albedo [20], calculation of radiative forcing[21], land 
cover classification, cloud detection, atmospheric 
corrections, and aerosol optical properties [18]. 

Accurate BRDF time series at customized spectra and 
spatial scales can estimate many biophysical phenomena 
that are currently wrought with errors. For example, up to 
90% of the errors in the computation of atmospheric 
radiative forcing, which is a key assessor of climate change, 
is attributed to the lack of good angular description of 
reflected solar flux and earth radiation budget (ERB)[22]. 
MODIS albedo retrievals show errors up to 15% due to its 
angular and spatial under-sampling when compared to CAR. 


Accuracy of BRDF estimation is therefore a representative 
metric of the ‘goodness’ of angular sampling. 

Angular Trade-offs with Spatio-Temporal Sampling 

Angular spread at the same ground spot affects spatial 
sampling which in turn affects temporal sampling. 
Increasing the Foresight angle of view elongates the sensor 
footprint and coarsens spatial resolution[6]. More the 
angular spread required, more will be the coarsening. This 
restricts the nadir resolution to a significantly small ground 
sample distance (GSD), to allow for resolution requirements 
for off-nadir pointing. Small GSD results in small swath, 
lower spatial coverage and less frequent revisits. Swath can 
be increased by adding more spatial pixels (and 
compromising on the spectral dimension [6]). Alternatively, 
swath can be kept the same and revisits increased by adding 
more clusters in the form of the clustellation. Each cluster 
would contribute to increasing angular spread while more 
clusters can improve spatial coverage and temporal revisits. 

4. Software Tools for Orbit Design 

Computational tools for spacecraft and system design have 
been very important in making early design decisions. 
Model-based systems engineering is a focus of working 
groups under the International Council on Systems 
Engineering (INCOSE) and the developed tools have been 
applied to existing missions such as RAX [23] and 
PHOENIX[4]. Tools for space logistics and interplanetary 
transportation such as SpaceNET, modular and open source, 
are also available [24]. The CubeSAT standard and 
associated documentation also provides a great resource to 
develop and integrate upto 3U (4 kg) spacecraft[25]. 

Existing Tools in Distributed Space Missions 
Existing tools for monolithic spacecraft and other space 
design can be and have been adapted for distributed space 
systems. Individual components of space system design can 
be combined from different software. For example, orbit 
design can be done using NASA GSFC’s GMAT (General 
Mission Analysis Tool) or NASA JSC’s Copernicus tool. 
Spacecraft operations can be aided by NASA JPL’s Activity 
Plan Generator (Automated Scheduling and Planning 
Environment (ASPEN) or Maestro tools. Specific interfaces 
for risk and science return for Saturn and Mars missions are 
also available. Tools for specific science data analysis such 
as USGS’s Integrated Software for Imagers and 
Spectrometers (ISIS) and ESA’s Rosetta Science Planning 
tool can be modified for some mission design. Cost/risk 
associated with distributed launches, staged deployment and 
reconfigurable constellations, all of which allow flexible 
design with increased costs, have been studied at MIT [26]. 
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Need for a new Tool for Earth Observation 
While all the above tools are great for specific missions and 
specific components, there is no off-the-shelf, modular tool 
for DSM design that can be used at the high-level 
architecture phase when key decisions are made. With the 
advent of hundreds of small satellites currently in orbit, and 
companies such as PlanetLabs and Skybox launching 
constellations in dozens, there is need for an integrated and 
modular tool which will enable easy plug into science 
metrics such as those for earth observation (extendable to 
astrophysics or navigation, etc.). Such a tool will allow rapid 
simulation of hundreds of architectures and their evaluation 
so that the “best” ones can be selected early in the design 
cycle. Since multi-angular earth observation has been 
established as an important field, we have developed such a 
tool for decisions related to DSM design for multiple angles. 


5. Relevant Formation and 
Constellation Models 

Formation flight has been analyzed at several levels of 
fidelity, focusing on understanding and manipulating the 
relative motion of satellites in the Local Vertical Local 
Horizontal (LVLH) frame, for achieving large angular 
spread for any spot under the cluster. Constellations are 
analyzed in the global. Earth-centered and fixed frame 
(ECEF) for achieving large angular spread for any/all points 
on Earth within acceptable time spans. 

Review of Formation Flight Models 
The linearized Hill, Clohessy and Wiltshire equations, 
simplified to be known as the Hill’s equations [27], [28] 
describe relative motion between any two spacecraft in a 
cluster, and can be extended to multiple spacecrafts. In this 
framework, one satellite is assumed to be traveling in a 
circular Keplerian orbit while the others are perturbed from 
this orbit by a small quantity. The HCW Equations [27], [28] 
has the X axis pointing radially away from the earth and Y 
axis in the direction of motion, is given by: 

= X — 3n^x — 2 ny 
ay = y 2 nx 

Uz = z -\- n^z 

Equation 1 

The additional orbit perturbations over and above these 
accelerations are J2 effects due to non-spherical Earth (~2.4 
X 10'^ m/s^ in LEO), third body perturbations due to 
differential force by the Sun and Moon on the spacecrafts 
(-3. 6-4. 3 X 10'^ m/s^ in LEO), solar radiation pressure 
(-1.7X 10'^^ m/s^ in LEO) and atmospheric drag due to 
small differences in the spacecraft shape and ballistic 


coefficient and atmospheric properties (-3.2X 10'^ m/s^ in 
LEO). The closed solutions to the Hill’s equations i.e. 
relative geometries which do not need any active control to 
keep them intact have an analytical form with 6 initial 
conditions. The formations of interest within these forms are 
the string of pearls (SOP) where the satellites remain in a 
string in the along-track direction separated by a constant 
distance, say S km, cross track scan (CTS) which is the SOP 
configuration extended to include oscillations in the Z 
direction of any amplitude and phase desired and the free 
orbit ellipse (FOE) where all the satellites arranged in 
elliptical rings around the LVLH origin. 


Since BRDF estimation requires inter- satellite zenith angles 
up to 80°, very large inter- satellite distances are required 
which violate the assumptions of the HCW equations. HCW 
does not account for Earth’s curvature and when propagated 
over a large time, the non-linear dynamics introduce large 
errors. Dual Spiral equations [29] provide relative equations 
of motion that analytically factor in the curvature of the 
Earth. They represent the motion of a point about a 
secondary axis which in turn rotates about a primary axis. 
P(Aa,5) represents the motion of the point in angles with 
respect to the primary axis, p 2 the angular radii of the point 
P from the secondary axis, (u2,52) the pole of the secondary 
axis with respect to the primary axis and co, the rotation rate 
of the secondary axis. 


S = 90 -a cos [sin S 2 cos + cos 5^ sin p^ cos cot\ 


= a cos 


2 


cos P 2 - sin 5^ sin 5 
cos 82 cos 5 


Equation 2 



Figure 2: [Right] Two satellites in co-altitude circular 
orbits separated by relative inclination (iR) and phase 
differences (c|)r). [Left] Relative motion of satellite 2 (S2) 
as seen from SI (orbital plane marked by arrow) in an 
earth centered coordinate system, full sky geometry. 


Under the assumption of the primary and secondary axis 
being mutually perpendicular at all times and co representing 
the rotation of the Earth, the dual spiral equations reduce to 
the relative analemma equations which describes large- 
scale relative motion of co-altitude satellites in circular 
orbits using two key parameters are important, relative 
inclination (iR) i.e. the angle at which the orbits of the two 
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satellites intersect, and the relative phase (cpR) i.e. the 
angular separation between them when one passes through 
the others’ place [29]. Relative inclination (iR) and relative 
phase (cpR) are: 

cos = cos cos I2 + sin sin I2 cos AN 
= (T2 — + AO 


AO = O2 + - Oi 

cos I2 - cos cos 

cos 0. = ^ ^ ^ ^ 

sin sin I2 

cos I 2 - cos cos 

cos 0 . = ^ ^ ^ ^ 

sin sin I2 

Equation 3 


Where AN is the angular separation along the equator of the 
two ascending nodes as shown in Figure 2a, T2- T1 is the 
time between the satellites crossing their respective 
ascending nodes, n is the angular motion of either satellite at 
their altitude and AO is an intermediate variable (Figure 2b) 
representing a difference in arc length from where the two 
orbits intersect to their respective ascending nodes. For 
circular co-altitude orbits, the relative motion of any satellite 
with respect to another is an analemma or a figure-of-8 
motion as seen in Figure 2b, left. Here, S2 is the base 
satellite and the analemma is the motion of SI as seen by 
S2. Parametric equations for the analemma in the inertial 
earth-centric frame are given by the following form where a 
is the azimuth of S2 at the earth’s center about the point 
where it crosses the orbital plane of SI and 6 is the elevation 
of S2 above the same point: 

sin S = sin sinnt 
a = nt — atan(cos iR tan nt) 

Equation 4 


In the same orbit plane then the analemma reduces to a point 
offset from the base satellite by an amount equal to the 
spacing between the two satellites. The nature and shape of 
the relative analemma trajectories indicate their usefulness 
for capturing both the SOP and CTS configuration by 
including the curvature of the Earth analytically. 

Higher fidelity models account for perturbations such as 
atmospheric drag, solar radiation pressure, non spherical 
earth and third body effects that accumulate over several 
orbits and need to be corrected for periodically. The 
modified HCW equations introduce the effects of J2 
perturbations due to the oblate shape of the Earth [30], and 
have an analytical form (Equation 5) with parameters 
described in [30]. When both J2 and atmospheric drag 
effects [31] are accounted for, the state matrix assumes a 
7X7 form (Equation 6), which is a modified version of the 
6X6 dynamic state matrix, and is required to be solved 


numerically to compute relative satellite trajectories. 
Parameters in Equation 6 are described in [31]. 
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For the highest level of fidelity, orbit modeling and 
propagation software such as Analytical Graphics Inc. 
Systems Tool Kit’s (AGI-STK[32]) High Precision Orbit 
Propagator (HPOP) is available. Full orbit propagation can 
be performed for all the satellites in the cluster and the 
solutions mapped into the LVLH frame with respect to a 
reference satellite. Closed cluster flight trajectories are 
formed when satellites are in orbits whose Keplerian 
elements (except the semi major axis) differ by a small 
amount or differentially. The relationships between the 
HCW coefficients in some simplified cases and the 
differential Keplerian elements is given by the COWPOKE 
equations or “Cluster Orbits with Perturbations of Keplerian 
Elements” [33]. The semi major axis is an exception 
because it corresponds to orbit energy, differing of which 
will break the formation. For circular LEO orbits, the HCW 
and differential Keplerian elements are related as: 

Ao = —aSe 

Bo = ay/ Si'^ sin^ iSCl^ 
x{t) = Aq cos{nt + a) a = 0 

y{t) = -2Ao sm{nt + a) + yoff f3 = uj - 9z 

z{t) = Bo cos{nt + (3) yof / = a{Suj A SM + cos iSQ) 

Equation 7 
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Review of Relevant Constellation Models 

The utility of constellations for multi-angular observation is 

two-fold: 

1. For global and/or more frequent coverage for NFOV 
payloads in formations as clustellations. 

2. For global coverage for WFOV payloads as 
constellations, which thereby provides angular coverage. 

Global, temporal sampling via constellations have been 
studied in the past [34] [35]. Only relevant literature to 
achieve multi-angle coverage will be focused here. WFOV 
constellations for MA-EO are most applicable to estimate 
the Earth’s radiation budget (ERB). Sun-synchronous orbits 
are not good for this application because they miss the 
extreme of the systematic diurnal variations. Walker 
constellations have been studied recently [34] [35] to be 
among the most efficient configurations for wide area 
continuous coverage. They rely on symmetric geometry and 
circular orbits to describe a constellation in only three 
variables corresponding to the number of planes, the number 
of satellites in each plane, and an inter-plane phasing 
parameter. Reference [36] and [37] has optimized Walker 
constellations to minimize global and regional revisit time 
respectively, and published the optimal number of planes, 
satellites and phasing for a given sensor field of view. 

Continuous and complete global coverage is also provided 
by the Streets of Coverage pattern [38], with excess 
coverage at the poles [35]. Flower constellations [39], [40] 
provide a more generalized framework for relative orbits 
and coverage than Walker constellations. The flower 
satellites have identical altitude, inclination and eccentricity, 
like Walker, and are characterized by them as well as three 
additional integer parameters. The authors prove that they 
can project any symmetrical shape onto the Earth Centered 
Earth Fixed (ECEF) just by varying the defining parameters 
[41], can make it rotating frame independent and make the 
configuration compatible with J2 invariant orbits because 
altitude -inclination are free design variables. They have 
been applied theoretically to Earth Observation [42], among 
others such as communication and navigation. 


6. Proposed Framework for Multi- 
Angular Tradespace Analysis 

The framework to assess the optimal formation architectures 
(unique combination of design variables such as orbit 
parameters, payload FOV, imaging mode, etc.) and validate 
their BRDF estimation capabilities couples Model-Based 
Systems Engineering (MBSE) with Observing System 
Simulation Experiments (OSSE)[ll]. A tradespace of 


formation architectures can be analyzed by varying the 
design variables in the MBSE model and assessing its effect 
on data assimilation and science products using OSSEs, as 
shown in Figure 3. Only the constellation and formation 
flight architecture generation and preliminary evaluation 
will be discussed in this paper. 



Figure 3: Summary of the overall approach to 
enumerate mission architectures under technical 
constraints and evaluate them based on proposed science 
products (right), tech specs and cost (left). 


^INPUT 

Orbits and 
Modes 
Module 
4 \ 


Payload 


Intersatelllte 

Comm 

Capabilities 


GNC 
Module 

^ ^0 reg Images Onboard 

for guidance 

Processing 

and Comm 


Number of sats 
Orbit Parameters 
Prob. Of Collision 


Position and 

Pointing 

Requirements 


INPUT 

BRDF Measurement Requirements 

• Spatial resolution (< 500 m) 

• Measurement Zenith Angles (< 60°) 

• Measurement Azimuth (< 360* *) 

• Solar zenith Angles (< 80°| 

• Spectral Resolution (>14 bands) 

• Spectral Range|3S0 2300 nm) 


Propulsion 

Module 


Nanosat Bus Requirements 

• Mass < 10 kg 

• Cube < 10X20X30 cm 

• Power < 25 W 
•Altitude (400- 800 km) 


OUTPUT: Metrics 


•Altitude 
OUTPUT: Metrics 
TRL, COTS support? 


OUTPUT: Metrics 


MPV, dV. 

Thrust, TRlt f 

Cost 

Module 


COTS support? 
Image D/L Rate 

OUTPUT: Metrics 
COTS support? 
Operation Life 


Figure 4: Systems Engineering Model (MBSE) 


The systems engineering model for multi-angular 
applications has been enumerated in Figure 4 as an N2 
diagram. Vertical arrows represent inputs into the modular 
subsystems (grey boxes) and horizontal arrows represent 
outputs. The model is iterative, as seen by the feedback 
loops and internal dependencies between subsystems. This 
paper concentrates only on orbits, or the first box. The 
outputs from the orbits module depend on if the DSM is a 
formation or constellation. The main output is angular 
sampling over time and temporal revisits (given an imaging 
mode for formations). The simplest imaging mode [12] for 
formation flight is assumed for this paper where a constant 
satellite is chosen as reference and looking nadir, while 
others point to the ground spot below the reference. Trades 
and design of the payload [6], subsystem capabilities to 
support such a mission[10], appropriate cost models [43] and 
coupling with appropriate OSSE models[2], [11], [12], [13] 
have been discussed in other literature. 
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Figure 5: BRDF data collected by the Cloud Absorption 
Radiometer during two NASA airborne campaigns. 


The science evaluation model is based on complex BRDF 
models and validated against data collected by the airborne 
Cloud Absorption Radiometer (CAR) instrument, developed 
and operated by NASA Goddard Space Flight Center 
(GSFC) [16] [20]. The CAR is designed to scan from zenith 
to nadir with a field of view (FOV) of 17.5 mrad and has 14 
bands between 335 and 2344 nm. By flying the instrument 
around a particular ground spot in circles and at different 
heights it is possible to get thousands of multi-angular and 
multi-spectral radiance measurements used for the accurate 
estimation of BRDF [44] [45]. Figure 5 shows the land 
cover type distribution globally, as extracted from the 
NASA MODIS database. For each type, the angular 
reflectance as measured by CAR is shown as a polar plot 
where radius is the measurement zenith and azimuth, the 
relative azimuth angle with respect to the sun. The 
reflectance shown is at representative wavelengths: 1032 nm 
for snow, 432 nm for water, 682 nm for grasslands and 
forests, 870 nm for croplands, cities. The local angular 
dependence of reflectance is apparent and the goal of our 
mission is to characterize this dependence globally, not just 
by an oversimplified biome type and from scarce air data. 

The output metrics mentioned in this paper will feed into the 
science evaluation model (described in [13] [2]) and be used 
to decide the final orbits for the formation/constellation. The 
method has been published for albedo products using 
formations [12] and ERB using constellations [1 1]. 

7. Implementing Distributed Orbits using 
THE Proposed Framework 

This section describes the modular implementation of the 
orbit generation module, both for formation flight and 
constellations. The architectures are enumerated to enable 


multi-angular earth observation so need large baselines (for 
NFOV angular spread) or multiple overlaps (for WFOV 
spread). The temporal trades are part of the WFOV angular 
trades because angular coverage per spot is achieved by 
satellite ground spot overlaps within an acceptable time. 



Figure 6: Process Flow Chart for tradespace analysis of 
DSMs with respect to constellation and formations. 


Figure 6 represents how constellations and their variables in 
the global, ECEF frame (left Box I) and/or formation flight 
geometries and their variables in the LVLH frame (right 
Box III) affect the angular and temporal performance 
metrics in the center (Box II). The formation variables are 
not listed because they depend on the models listed (will be 
detailed in the next section). The metrics in Box II feed into 
the science evaluation or OSSE model. There are some 
formation imaging modes which use the constellation 
analysis software because they are optimized using 
performance at the ECEF ground spots (not in LVLH). The 
software’s versatility to enable this demonstrated, but the 
application described elsewhere [12]. 

Formation flight 

Formation flight for multi-angular measurements (Figure 1- 
left) can be analyzed in a 3 -layer framework of increasing 
model fidelity and decreasing computational ease of 
exploration - Figure 7. Each level’s trades are analyzed in 
detail to streamline and inform the selection of variables in 
the next (higher fidelity) level. The final metrics are 
generated from M ATE AB -controlled STK, however the 
streamlining exercise helps control the explosion of design 
variable combinations in a computationally expensive 
environment. 



Figure 7: Levels of models used for formation flight 
simulation as a function of fidelity and computational 
ease of tradespace exploration. 
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All analysis in this section has been discussed in the Local 
Vertical Local Horizontal (LVLH) frame centered at the 
reference satellite at a 500 km altitude. The first level is the 
linearized Hill Clohessy Wiltshire Equations (Equation 1) 
which are numerically corrected to account for Earth’s 
curvature at very large inter satellite distances by 
transforming from Cartesian to curvilinear coordinates. The 
HCW solution that gave the most angular diversity at the 
ground target was the Free Orbit Ellipse (FOE) 
configuration demonstrated in Figure 8 with 12 satellites + 1 
looking nadir. The tradespace variables are the ring radius of 
the ellipses, their shape and their inclination with respect to 
the chief orbit. 



Figure 8: Free Orbit Ellipse simulated using the Hill 
Clohessy Wiltshire (HCW) Equation 1 with 12+1 
satellites and numerically corrected for Earth curvature. 

The satellite trajectories are in LVLH and the orange 
star represents the LVLH reference nadir. 

Figure 9 shows the trajectories (top) and the angles they 
subtend on the ground (bottom) when the initial x and z 
positions of the satellites are varied to give different 
orientations to the chief orbit. Relative azimuth and solar 
zenith angles are not plotted for simplicity. The red and 
green curves correspond to ellipses that project circles on 
the ground and are circles respectively. Since the trajectories 
represent the movement of each satellite over on orbit and 
the biome/ground spot below the cluster changes 
throughout, the mission designer can choose the optimal 
orientation based on the type of angular spread per biome 
they want. The COWPOKE equations say this implies 
tweaking the differential eccentricity, RAAN (Equation 7). 

The cross-track scan (CTS) and string of pearls (SOP) 
solutions were achieved using the dual spiral equations 
(Equation 2) simplified into the relative analemma equations 
(Equation 3). These solutions analytically accounts for the 
earth’s curvature, azimuthal variation is not as easy to 
achieve as the FOE solution. The analemma equations can 
be customized for BRDF -related relative motion by 
transforming the parametric equations (Equation 4) from 
inertial earth centric coordinates to LVLH coordinates using 
the knowledge that the chief orbit is at a distance of (Rn+h) 


from the center of the earth, h being the orbital altitude and 
Re, the radius of the earth. The transformation is given by: 

(n 


Xk(t) = i^E + h) 
yfc(t) = (Re + h) 

Zfc(t) = (Re + h) 


cos(Ofc + a(t)) sin 


5(t) 


sin(0^. + a(t)) sin 
n 


m 


cos 


m 



Equation 8 
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Figure 9: [Top] Curvature corrected FOE trajectories in 
LVLH frame for changing ellipse inclinations to chief 
orbit (xO/zO ratio); [Bottom] View zenith angle 
subtended by a satellite at the ground spot directly below 
the LVLH origin for the FOE configurations shown. 

Equation 8 represents the motion of the k’th satellite, 
located at a phase separation of (= in Figure 2) in 
the LVLH/HCW frame as seen by a base satellite that may 
be real or virtual, always located at the origin of the HCW 
frame. A 9 satellite and 1 reference satellite case in the CTS 
configuration (axes not equal) is shown in Figure 10- 
trajectories in blue and their projections on the 
perpendicular planes in red. The differential inclination and 
phases may be varied to get a large option of angular 
spreads over one orbit, as seen in Figure 11. The three 
curves toward the bottom correspond to those with 
negligible differential page (middle trajectories in Figure 
10). The six curves at the top of Figure 11 are the flanking 
analemmas, the more the differential inclination, more the 
size of the analemma. Maximum azimuthal coverage is 
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obtained when the satellites go toward the extremes of the 
analemmas twice every orbit. Each satellite will have a 
unique analemma without any overlapping phase with 
another, to prevent collisions at the highest latitude, unlinke 
FOE where multiple satellites can share an ellipse. 
Verifying with the COWPOKE equations (Equation 7), for 
circular orbits, the only variable to vary azimuthal spread at 
a given phase is differential inclination. 

The SOP configuration provides the same angular coverage 
as a monolithic satellite with many forward aft sensors (e.g. 
MISR) and provides limited azimuthal coverage, hence has 
not been explored in much detail. It can be simulated using 
the analemma equations but with no differential inclination. 



Figure 10: Cross Track Scan (CTS) simulated with 9+1 
satellites with differential inclination and phase using the 
relative analemma model (Equation 8). The satellite 
trajectories are in LVLH and the orange star represents 
the LVLH reference nadir. 



Figure 11: Variation of the view zenith angle subtended 
by a satellite at the ground target directly below the 
LVLH origin for the 9 satellites in Figure 10. 

Since the HCW and dual spiral level of analysis does not 
account for perturbations, the trajectories for every orbit are 
exactly the same. The second level of analysis introduces 
the effects of J2 perturbations due to the oblate shape of the 
Earth [30](Equation 5) and then atmospheric drag effects 
[31]. Figure 12 shows the trajectories of 3 satellites in 
different colors with different ring sizes, simulated by 
maximum X-intercepts (Equation 6). The orange star is the 


ground target directly below the reference satellite (marked). 
Initial X and Z positions and velocities can be varied in 
keeping with the HCW equations to get a large tradespace of 
trajectories corresponding to the HCW variables of ring 
sizes, shapes and orientations. 


Initial X = 100 km 

Initial X= 160 km 



Figure 12: Free Orbit Ellipse with 3+1 satellites, 
simulated using the modiHed HCW equations with J2 
and drag (Equation 6 with appropriate initial 
conditions), propagated for 1 day. The satellite 
trajectories are in LVLH and the orange star represents 
the LVLH reference nadir. 

The expected tumbling effect of free-orbit ellipses about the 
cross track axis due to J2 effects is clearly seen in Figure 12. 
Additionally, there is a slight drift in the along-track 
direction due to atmospheric drag (negligible at 500 km 
altitude). The corresponding view zenith angles for all three 
rings, simulated at three different orientations, are seen in 
Figure 13. More the initial X position, larger the ellipse and 
more the angle subtended at nadir, as expected. The ellipses 
that make a 45 deg or higher angle with the LVLH horizontal 
display two crests per orbit. The crests correspond to the 
two extreme positions in the along-track direction. They are 
symmetric at the beginning and lopsided after a day due to 
J2 tumbling in that direction. The zenith angle minima 
corresponds to the higher altitude side of the ellipse and the 
trough between the crests the lower altitude side. As the 
initial X-Z position ratios are varied to increase the 
inclination angle of the ellipse, the double crests become 
single because the ellipse tends toward a horizontal circle. 

In terms of the angular output, lower inclinations are useful 
when a large variation of view zenith angle is needed over 
the orbit, for example, to follow the variation of the solar 
zenith angle and precisely estimation specular reflection or 
hotspots. Higher inclinations are useful when an 
approximately constant view zenith is desired. Note that the 
relative azimuth with respect to the sun will show the same 
variation for all ellipses, because they share the same plane. 
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Figure 13: View zenith angle subtended at the ground target directly below the LVLH origin for the 3 satellite 
trajectories in Figure 12 (corresponding colors). The different line types correspond to different inclinations of the 
FOE with the LVLH horizontal plane containing the chief orbit (only one in Figure 12) 


To introduce azimuthal variation, an eccentric chief orbit is 
required with differential Keplerian[33]. The corresponding 
equations (also COWPOKE) are more complex than 
Equation 7 and have not been analyzed in this paper. 
Parallel literature [10], [12] has shown that slightly eccentric 
orbits are beneficial in compensating for the drift in the 
along track direction and aid in keeping the formation 
together after 7-8 months of operation. It would be valuable 
to find the optimal chief and differential eccentricity that 
would allow a good azimuthal spread via free orbit ellipses 
as well as be beneficial to maintain. This analysis has been 
assigned to future work, alongside determining optimal 
chief orbits for the formations discussed in this paper (sun 
synchronous vs. frozen vs. polar vs. eccentric). 

The 7X7 state transition matrix can be initialized with a 
finite Y and Z position, in order to simulate cross track 
scans that account for J2 and drag effects. As before, a 
constant atmospheric density of 0.02 kg/m^, and a 6U 
cubesat form factor for the satellites was assumed. Figure 14 
shows 3+1 simulated trajectories, propagated for 1 day, 
where the reference satellite is at (0,0,0), directly above the 
orange star. The pattern corresponds exactly to Figure 10, 
except plotted on equal axes to differential between the 
extent of radial and along track drift due to drag and J2. The 
latter is obviously far more, however it does not affect the 
formation because the relative spatial and angular spread 
between the satellites does not change much in a day. 
Global simulations using STK will show how these drifts 
are hugely affected by the differential Keplerian elements 
used to achieve the initial conditions. 



Figure 14: Cross track scan with 3+1 satellites, 
simulated using the modified HCW equations with J2 
and drag (Equation 6 with appropriate initial 
conditions), propagated for 1 day. The satellite 
trajectories are in LVLH and the orange star represents 
the LVLH reference nadir 

The third and last framework of models uses AGFs Satellite 
Tool Kit to initialize and propagate individual satellite orbits 
(High Precision Orbital Propagator or HPOP) and then 
calculates their relative trajectories with respect to a 
reference satellite. SOP, CTS and FOE configurations of 
varying shapes, sizes and orientations can be created by 
varying the differential Keplerian elements of the satellite 
orbits. For example, verified by COWPOKE[33], when 
multiple satellites have the same Keplerian elements except 
separated by a small, differential true anomaly or TA, the 
resultant relative motion is the string of pearls (SOP). When 
multiple satellites have the same Keplerian elements except 
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separated by a small, differential TA as well as inclination 
or RAAN, the resultant relative motion is the cross track 
scan (CTS). Finally, if satellites have different TA, 
inclination or RAAN as well as perigee or eccentricity, the 
resultant relative motion is the free orbit ellipse (FOE). 
External constraints like biome and latitude of interest 
determine the range of some variables for the angular trades. 
For example, differential inclination (RAAN) produces 
maximum separation at the poles (equator). FOE and CTS 
should be created by either one depending on target latitude 
where maximum angular spread is desired. 

Figure 15 shows the trajectories of 3+1 satellites with 
differential inclinations and TA, propagated on STK for 1 
day. The drift associated with the formation to achieve 
hundreds of kms of baseline and angular spread (like Figure 
10) is significant. Different inclinations cause the J2 forces 
on each satellite to be different causing the RAAN to rotate 
differently. A 20 deg view zenith angle spread takes less 
than 3 months to break up [10] because increasing RAAN 
differential causes the angular spread at the poles to 
decrease and equator to increase and eventually the latter is 
too large for the satellites to see each other. Figure 15 
shows a drift not only in the along-track direction (as seen in 
Figure 14) but also in the radial direction, indicative of the 
formation breaking (similar to semi major axis difference). 



Relative In-Track Motion 


Figure 15: Cross track scan with 3+1 satellites, simulated 
on STK using differential TA and inclination, 
propagated for 1 day and satellite trajectories exported 
and plotted in the LVLH frame. The orange star 
represents the LVLH reference nadir. 

When free orbit ellipses is simulated on STK using 
differential inclination, with TA and eccentricity, 
propagated for a day and plotted in the LVLH frame, a 
significant along track drift and tumbling effect is seen - 
Figure 16. This observation can be attributed to increasing 
RAAN spread due to the very different inclinations required 
to achieve the large baselines. The full extent of this drift is 
not seen in the modified HCW level of analysis. The loss of 


symmetry in the two-crests-per-orbit pattern for VZA is 
much faster in this model (Figure 17) than Figure 13. This is 
because there is no way to distinguish between LVLH-Z 
(cross track) due to inclination vs. RAAN. In reality, both 
have significant contribution, as is apparent from the 
COWPOKE equation for Bq (Equation 7), and can be used 
for design. The modified HCW analysis is thus used as an 
informative step in streamlining the maximum prescribed 
baseline for maintainable clusters, but the exact numbers 
and the ways to achieve those baselines are not relied upon. 


M4t0n tjf Im wi1 a i = TA : 0 i) = 0 



Figure 16: Free orbit ellipse with 3+1 satellites, 
simulated on STK using differential RAAN, TA and 
inclination, propagated for 1 day and satellite 
trajectories exported and plotted in the LVLH frame. 
The orange star represents the LVLH reference nadir. 
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Figure 17: View zenith angle subtended at the ground 
target directly below the LVLH origin for the 3+1 
satellite trajectories in Figure 16. 


As indicated by the COWPOKE equation for Bq (Equation 
7), cross track scans can be simulated using differential 
RAAN and TA. 3+1 satellites were propagated for 1 day on 
STK and the relative trajectories plotted in Figure 18. It is 
immediately apparent that compared to Figure 15, the drift 
in all directions is far less even with hundreds of kms of 
baseline and hence great angular spread on the ground. 
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When the free orbit ellipse is simulated using differential 
RAAN instead of inclination, similar improvement in drift is 
seen (Figure 19). Some tumbling effect is seen as expected, 
in keeping with Figure 12. Going forward, only differential 
RAAN and TA (up to 6 deg) will be used to achieve 
formation spreads. 6 deg shows manageable drift, and 
sufficient angular spread, as will be shown in the next 
section. Since such formations will provide no azimuthal 
spread at the poles, different imaging modes (strategies to 
point the payload) will be used to achieve angular coverage 
at the poles, as described in other literature [12]. 



Figure 18: Cross track scan with 3+1 satellites, simulated 
on STK using differential RAAN and TA, propagated 
for 1 day and satellite trajectories exported and plotted 
in the LVLH frame. The orange star represents the 
LVLH reference nadir. 



Figure 19: Free orbit ellipse with 3+1 satellites, 
simulated on STK using differential RAAN, TA and 
eccentricity, propagated for 1 day and satellite 
trajectories exported and plotted in the LVLH frame. 
The orange star represents the LVLH reference nadir. 

BRDF has three major angles of interest, as mentioned 
before. Only view zenith angles have been shown in the 
paper so far. Solar zenith angles for the cluster and the 
relative azimuth angle with respect to the sun for every 
satellite in the cluster have been analyzed but not presented 
for brevity. It is important to note that those angles are 
equally important, and all three angles, i.e. 2N+1 angles for 
an N-satellite cluster per instant of time, serve as outputs 


from the orbits module and inputs into the science 
evaluation module driven by OSSEs. 

To demonstrate angular spread in the global frame, 
three candidate clusters with 9 satellites each (to match 
MISR’s sensor numbers) were simulated in STK to image a 
specific spot on earth [0, -103.729] at a repeat period of 16 
days and compared to the measurement spread of the same 
ground spot by MISR. The arrangement is seen in Figure 
20-top. One of the three clusters was in the SOP 
configuration (black) where in all Keplerian elements except 
the true anomaly were the same. The two other clusters were 
in different FOE configurations (blue and green), simulated 
by perturbing differential RAAN, eccentricity and TA only. 
There was an approximate 60° RAAN differential between 
the chief orbit of MISR and any of the clusters. Due to 
MISR’s large swath, it has a total of ~ 1800s of access to the 
ground spot while the clusters have only 1.4s in 16 days. 



Zenith in <fegr?s^ 

Figure 20: [Top] Global simulation of clusters (green 
dots as different architectures) vs. the MISR instrument 
with its 9 sensors (pink) [Bottom] BRDF polar plot for 
simulated measurements made by three 9-satellite 
clusters in SOP (black) and 2 different FOE 
configurations (green, blue) and the MISR instrument 
using 9 sensors (pink) of the same geographic location 
shown in (a). STK was used to calculate the 
instantaneous view zenith (radius) and view relative 
azimuth (azimuth) angles. The solar zenith angle is ~89° 

The simulated angular measurements of the target are 
plotted on a BRDF polar plot for the time instant when the 
reference satellite is directly overhead the target - in Figure 
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20-bottom. The radius shows the view zenith angle and the 
polar azimuth the relative azimuth angle with respect to the 
sun, for a constant solar zenith. Since MISR has fixed 
sensors, the measurement zenith angles for a direct overpass 
(pink asterisks in Figure 20-bottom) are pre -determined and 
the relative solar azimuth is the angle between the velocity 
and sun vector measured in the satellite HCW X=0 plane. 
These measurements for a single overpass can be easily 
replicated by an SOP cluster (black asterisks in Figure 20- 
bottom). Improvement in angular performance is clearly 
demonstrated in both the FOE clusters because a large and 
diverse azimuthal and zenith spread on the BRDF plot is 
possible. Different architectures of different combinations of 
number of satellites and their differential Keplerian 
elements output such a spread at every instant of time, and 
these angular-temporal spreads serve as inputs to be 
evaluated in terms of how well they capture the ‘true’ BRDF 
in Figure 5[2], [12]. 

MATLAB -based software has been developed to automate 
formation architecture generation on AGI’s STK, angular 
report creation for each architecture and angular metric 
calculation as a post-processing step, with the option to 
evaluate as is or pass into a science evaluation model. 
Architectures are generated by permuting allowable 
combinations of design variables - number of satellites, 
chief orbit altitude, inclination and differential RAAN and 
TA and the reference satellite among them. To prevent 
variable space explosion and in keeping with the previous 
streamlining analysis, the number of satellites is limited 
between 3 and 8. Three satellites is the minimum required 
for the BRDF OSSE models and eight corresponds to 
NASA arc’s Edison Demonstration, currently the highest 
number of commissioned satellites in any DSM. For any 
given number of satellites (say, N), N-1 RAAN-TA 
differential combinations are picked from the 8 available in 
Table 1. These 8 slots have been selected as the corners of a 
+/-5 deg square in RAAN and TA, with the reference 
satellite in the center and no common TAs to avoid 
collisions. 


RAAN 

0 

-5 

5 

0 

-5 

5 

5 

-5 

TA 

-5 

-6 

-4 

5 

6 

4 

-1 

1 


Table 1: RAAN-TA (in deg) slots for the full factorial 
enumeration of formation architectures 


For a given altitude-inclination combination, that gives 1254 
RAAN-TA combinations for 3-8 satellites. Figure 21 
captures the software process flow. MATLAB -driven STK 
generates customized reports (some listed in the middle 
column of Figure 21) which are named in keeping with the 
design variables for easy post-processing. The ones in red 
boxes contain the measurement zenith and relative azimuth 
with respect to the satellite’s momentum vector, correspond 
to a 650 km altitude, 4 satellites, RAAN-TA combination #1 
and the reference sat number (the one that looks nadir). 
Each report (.csv) contains view (zenith or azimuth or solar 
with respect to 2 perpendicular axes) angles of every 
satellite with respect to the reference. These reports are post- 
processed to output the 3 angles of interest (VZA, SZA, 
RAZ) at every instant of time in the LVLH frame. 

Constellations 

When the DSM being analyzed is constellations, the 
analytical framework illustrated on the left side of Figure 6 
is used. Spatial coverage and sampling in the global ECEF 
frame, not the LVLH frame, are orbital outputs when given 
inputs from the payload module (for pixel sizes and 
spectrometer type). Parameters such as the grid size on the 
earth (default: 5degX5deg in both latitude and longitude) 
and time sampling (default: 1 minute) can be defined. Using 
these requirements, automatic scripts on MATLAB drive 
STK to generate multiple architectures on STK by 
permuting the orbit design variables. For example, three 
architectures are pictured in Figure 22 ’s left column. 
MATLAB -driven STK then commands each architecture 
definition to generate a full access report as a .cvaa file, 
some of which are seen in Figure 22 ’s second column. The 
reports are named as before, however with different design 
variables. For example, the first one is a 1 plane, 1 satellite 
design at a 400 km altitude/60 deg inclination, a 15 deg 
instrument field of view for latitude coverage below 40 deg. 



^ measAngles_alt650_nsats4_RTl_ref2.csv 
^ mcasAngles_alt650_n5ats4_RTl_reG.csv 
^5j^measAngl^^ltf5^jsats^T^ef4xs^_^ 
I ^ momenAngles_alt650_nsats4_RTl_refl.csv | 
momenAngles_alt650_nsats4_RTl_ref2.csv 
^ momenAngles_alt650_nsats4_RTl_ref3.csv 
^ momenAngles_alt650_nsats4_RTl_ref4.csv 
^ solAngles_alt650_nsats4_RTl_refl .CSV 
^ solAngles_alt650_nsats4_RTl_ref2.csv 
^ solAngles_alt650_n5ats4_RTl_ref3.csv 
^ solAngles_alt650_nsats4_RTl_rcf4 .CSV 
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Figure 21: Process Flow for Angular Metric calculation in the LVLH frame using different formation flight 

architectures, pre-defined by design variables 


13 





Lat40.fovl5.AltW0Jr»c«).Pldf»e$l.S*tsl.cvM 


Alt400_lr»c60_Ptanesl. 

Alt400.1nc60.Planesl. 

Alt400.Ii>c60.Pl*n«$l. 

AltW)Jr»c60_PUi>etl. 

Alt400Jr»c60_PUr»esl. 

Alt400>c60.PUnesl. 

Alt400.Inc60.PUnesl. 

AK400Jr»c60_Pldned. 

Alt400>c60_Plan«sl. 

Alt400>c60_Plan«sl. 

AMnn PI 


.Sa^.cvM 

.SatsS.cvsa 

,S«ts7.cvM 

.S«ts9.CVM 

,S«tsl.cvsa 

,Sats3.cvaa 

,S«ts5.cvM 

,S«ts7.cv»a 

,Sats9.cvaa 

.Satsl.cvaa 


Alt400Jr»c60_Pl«n«l_S«ti5.cvM 


Alt400_lnc60_Planesl. 

Alt400_Inc60_Plan€sl. 

Alt400.Ii>c60.PUn«sl. 

Alt«0Jr»c60_PUn«sl. 

Alt400Jr>c60_Plar»esl. 

Alt400>c60_PLin€sl. 

AhWJncfiO.PUnol. 

Alt400Jr»c60_Pldne$l. 

Alt400_lr»c60_Pl8nesl. 

Alt400>c60_Planesl. 

Alt400Ji>c60.PUn€il. 

Alt400Jr»c60_Pldne$l. 


,Sats7.cvaa 

.Sat^.cvaa 

.Satsl.cvM 

.SatsS.cvM 

,S«ts5.cvsa 

.Sats7.cvaa 

,Sitj9.cv*a 

,S«tsl.cvaa 

,Sats3.cvaa 

.SatsS.cvaa 

,Sdt$7.cvM 

,Sats9.cvaa 


/ 



Figure 22: Process Flow for Temporal Metric calculation globally or specific areas or latitudes of interest (ECEF 


frame) using different constellation architectures, pre-defined by design variables. 


Each access report is a detailed text file - Figure 22 third 
column - containing the time period (from when to when) at 
which every grid point is accessed by every sensor. These 
.cvaa files (per architecture) are post-processed to provide 
customized temporal metrics such as revisit time, time for 
first access, number of accesses, time for global coverage, 
etc. for every grid point and the globe. Science metrics 
dependent on such temporal metrics can easily be calculated 
by plugging in the science evaluation module, modeled after 
an OSSE. Since the analysis is global and across any length 
of mission lifetime, the metrics calculated are spatial and 
temporal. 

A separate script is also available for analyzing 
specific grid points in the same way as above, thus saving 
the computational resources required for full global analysis. 
For any target location or ground station, a text access report 
is automatically generated by M ATE AB -commanded STK 
and automatically repeated for multiple architectures. Each 
generates an access report, which can be post-processed to 
generate customized metrics as described above. 

STK, without the parallel processing license and dozens of 
available cores, was found to be inefficient for architecture 
studies using customized angular metrics. For example, the 
following angles are required to be calculated at every grid 
point and time instant for every architecture: MZA for each 
satellite (angle between the satellite vector at the ground 
spot and the zenith; <90 deg), SZA (angle between the sun 
vector at the ground spot and the zenith; <90 deg for solar 
spectrum) and RAZ (angle on the horizontal between the 
satellite vector projection and sun vector projection at the 
ground spot; <360 deg). The number of calculated angles is 
the product of the number of ground points (1651 by 


default), number of satellites, number of angles (3 by 
default), number of architectures and number of time steps. 
STK calculates all of them using a tool called Grid Inspector 
which loops over the number of points, satellites and angles, 
and re-calculates access for all loops, therefore taking 5 days 
to compute only 1 angle for a 64 satellite constellation. 

To improve efficiency, our tool uses STK only for temporal 
and spatial analysis, and contains a standalone module for 
customized metric calculations such as angular outputs. The 
algorithm is as follows: (1) The High Precision Orbit 
Propagator (HPOP) using the Jacchia-Roberts Atmospheric 
Model with up to J4 terms is used to propagate the all 
satellite orbits in every architecture and the resultant states 
per time step saved as text files. (2) The access reports for 
all architectures - as seen in Figure 22 - are saved for exact 
global coverage. (3) The grid point information is saved. 
The above three outputs from STK are then post-processed 
to calculate the required hundreds of thousands of angles 
offline. The results are validated against those calculated by 
STK’s Grid Inspector for one satellite propagated over one 
day. A reasonably good fit is seen for all grid points and all 
times with less than 5 deg. of average error. This error is 
less than half the angular resolution available by a grid and 
time resolution of 5 deg. and 1 min, for a satellite ground 
velocity of 7.3 km/s, hence considered negligible. All 
sampled angles, their dependent metrics or any other 
customized metrics dependent on global or temporal 
coverage, will be calculated in the above way. The orbits in 
our study are limited to LEO (the tool is unrestrictive). All 
other variables are set based on the case study. For example, 
in an ERB application, an FOV of 130 deg and a maximum 
of 64 satellites were simulated[ll]. 
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Figure 23: Validation of the VZA, SZA, RAZ calculated by the proposed tool (blue) at every grid point and time 
instant those calculated by STK-GI (red). The top panel shows one angle as seen by a single satellite at every minute 
for a day while the bottom panels show all three angles over a 22 minute period, for better visualization. Average 
error <5 deg is negligible for a grid resolution of 5 deg, time resolution of 1 min, satellite ground velocity of 7.3 kms/. 


8. Selected Results from Case Studies 

The analysis framework described is used to streamline the 
design variables, followed by the formation or constellation 
architecture generation software used to generate thousands 
of relevant designs for multi-angular observation. This 
section shows the angular and temporal outputs for some 
architectures, to demonstrate the utility of the tool in making 
decisions after considering conflicting objectives of 
performance. The objectives can be extended to spatial, and 
spectral sampling, and also cost. 

Formation Flight 

A 4-satellite formation with a circular chief orbit at 650 km, 
51.6deg inclination is considered. If the 3 non-chief 
satellites were to have one of the differential RAAN-TA 
slots in Table 1, 56 architectures are possible. The angular 
outputs (VZA, RAZ only) of two of those architectures are 
plotted in Figure 24 and Figure 25. The TAAN-TA 
combinations are mentioned in the legend. A 4-hour 
simulation in LVLH allows for >2 orbits and all-biome 
sampling. For formation maintenance analysis, >1 year 
simulations are required. The black lines correspond to 


MISR’s 9 cameras. The MISR instrument is onboard the 
TERRA monolithic spacecraft. Since it has 4 forward and 
back cameras at the same boresight angle, there are four, 
unique and constant VZA and two RAZ. 

For any relative orbit architecture, any of the satellites could 
act as a reference and point straight down while the others 
point to the spot below it. The colors in Figure 24/25 
represent the VZA and RAZ of each satellite in the cluster 
when a different one acts as reference. Figure 24/2 5 -top 
does not show 4 VZA curves per color because of overlaps 
due to symmetry in differential Keplerian elements. For 
example, the VZA curve for sat#3 with sat#l as reference 
(one of the blue curves) will be the same as the VZA curve 
for sat#l with sat#3 as reference (one of the cyan curves). 
Moreover, if the differential RAAN and TA of 2 satellites 
are symmetric, for example Sat#3 and #4 in Figure 25, it 
results in more overlaps. The solar zenith (SZA) for all 
architectures is nearly the same, because the satellites are 
fairly close and traverse nearly the same ground track with 
respect to the sun in the sky. The solar azimuth is very 
different, which contributes to the very complex RAZ plots. 


View zenith angle of each sat in the 4-sat cluster 



Figure 24: View zenith angle (top) and relative azimuth 
wrt Sun (bottom) of every satellite in a 4-satellite cluster 
in configured as one architecture - legend - when each is 
used as reference (grouped by color) and MISR (black). 

Bottom panels show the VZA-RAZ spread on a polar 
plot for t=50min for MISR (black) and cluster with 
Sat#l (blue) and #3 (green) as reference. 

Obviously, it is quantitatively impossible to judge the better 
architecture among the Figure 24/25 (or thousands more) 
looking at (2N+1) angle curves per architecture and N-times 
more if different reference satellites are considered. While 
polar plots (Figure 24-bottom panels) offer better 


representation of the angular spread at any instant of time, 
each N-sat architecture has N of these spreads for every time 
instant and solar zenith. The panels also show how different 
the spread is when a different reference is used for the same 
architecture. The only way to evaluate these angular spreads 
is thus to feed them into an OSSE which calculates the 
BRDF-error (or error in relevant products like albedo) per 
time step. The coupled MBSE+OSSE model is thus an 
indispensible judge to distinguish architectures based on 
science performance of their angular sampling. 


While NFOV formation angular outputs can be analyzed in 
the LVLH frame above, the WFOV angular outputs are best 
analyzed in the ECEF frame using constellation software. 



Minutes of simulation on October 11 2013 


Relative Azimuth angle of each sat in the 4-sat cluster 



Figure 25: View zenith angle (top) and relative azimuth 
wrt Sun (bottom) of every satellite in a 4-satellite cluster 
configured in another architecture when each is used as 
reference (grouped by color) and MISR (black). 
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Constellations 

For temporal output demonstrations, a Walker constellation 
with varying number of satellites and FOVs were simulated 
for a constant altitude of 709 km and inclination of 98.18 
deg., in keeping with the orbit of the A-Train and EOS 
satellites. Plane-Sat arrangement was not considered 
because it does not affect revisit time as long as uniformly 
arranged. Only latitudes below 70 deg. were considered. 
Figure 26 shows the maximum revisit time provided by all 
the architectures. A monolithic spacecraft with 15 deg FOV 
(e.g. Landsat) provides a 350 hour revisit - black circle on 
Figure 26 ’s left and full global map on Figure 26 ’s right top. 
The results show that at least 14 such satellites are required 
for a daily revisit (Figure 26 black line) and 16 satellites for 
a daily repeat (analytical calculation). Doubling the FOV to 
30 deg allows the same revisit in about quarter the satellites 
(4 satellites). The right panel also shows that revisits are far 
more frequent at higher latitudes than lower ones for polar 
constellations. 

While revisit time does not depend on constellation 
arrangement, metrics such as time required for full global 
coverage does. Figure 27 shows the time taken for the last 
grid point on the globe is accessed by different constellation 
architectures. The results from our tool show that global 
coverage is faster if the same number of satellites is 


arranged in more planes. In fact, lesser number of satellites 
(e.g. 8 satellites in red vs. 12 satellites in grey) can achieve 
coverage faster if arranged in more planes. The trade-off 
however is in terms of cost because launching into 8 planes 
requires 7 plane changes (each costing a bulk of fuel) or 
needs 8 times the number of launches than launching into 1 
plane. The increased performance and cost saved in 
developing 4 extra satellites can be compared against the 
cost of launching into 8 instead of 1 plane for the optimal 
design decision. 

The 8 satellite arrangement can be further analyzed by 
plotting the timeline of global access (until 100% is 
accessed) - Figure 28(bottom). Walker Delta arrangements 
are seen to be better than Walker Star in time to global 
access and all curves are compared to the monolithic 
counterpart in the same orbit. The monolith takes 14 days 
for full global access as seen earlier and is shown with a 
black line, called ‘Landsat’ because it is an existing 
spacecraft in the same orbit with a 16-day repeat cycle. 
Figure 28 ’s curves are very useful in deciding which design 
to choose depending on the integration time available over 
measurements and the coverage flexibility. For example, if a 
2 day integration time is available, then the 8 plane 
constellation is equivalent to having continuous global 
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Figure 26: Results from our temporal trade tool using Revisit time as a metric. [Left] Maximum revisit time over all 
gridpoints as a function of payload FOV and number of satellites in a Walker constellation at 709 km altitude, 98.18 
deg inclination. The thick black line indicates at least a 24 hour revisit for any point and the black circles indicate the 
designs for which global revisit time is shown on the right. [Right] Average revisit time at every grid point, calculated 
over a 16 day period, for a 1 (top) and 16 (bottom) satellite constellation. 
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coverage and there is no value in adding more satellites. 
Better angular output using WFOV sensor constellations 
and its utility in reducing uncertainties in the Earth’s 
Radiation Budget has also been demonstrated using this 
tool[ll]. It can thus be used to select a design, given a 
required temporal resolution, spatial or angular coverage. 


16 pUnM, 1 satellite 
8 planes, 2 satellite 
4 planes. 4 satellites 
2 planes. 8 satellites 

1 plane. 16 satellites 
12 planes. 1 satellite 

6 planes. 2 satellite 
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Figure 27: Time required to global coverage for Walker 
constellations (at 709 km, 98.18 deg) with varying 
number of satellites and their planar arrangement. 
Walker Delta constellations, on an average, show lesser 
time than Walker Star. 



Figure 28: Percentage of the globe covered with respect 
to time for Walker constellations (at 709 km, 98.18 deg) 
with 8 satellites arranged in different planes (colors) and 
in Delta (continuous line) or Star (dotted line) 
arrangements. 


9. Summary and Future Work 

This paper identifies a lack in open-source, orbital analysis 
software for understanding the complex design trades for 
distributed space missions (DSM), both constellations and 
more pertinently formation flight. A gap in the measurement 
making abilities of monolithic spacecrafts in the field of 
multi-angular earth observation is also identified and DSMs 
proposed to fill this gap. BRDF and its dependent products, 
such as albedo and ERB, are assumed to be appropriate 
metrics for angular sampling on Earth (solar spectrum). 

A framework and relevant software has been developed to 
analyze the full tradespace of LEO formations (in the LVLH 
frame over all Earth biomes) and constellations (in the 
ECEF frame) so that NFOV and WFOV sensors, 
respectively, which can provide angular coverage of all/any 
point on Earth. All software is written on MATLAB, STK 
or Excel and is completely automated in terms of 
architecture generation and evaluation. The formation 
models are arranged in increasing order of fidelity and 
computational requirements, and serve as a method to 
streamline the design variables as the models get more 
complex. The highest fidelity models are coded in 
MATLAB -driven STK and can generate thousands of 
architectures permuted from the streamlined variables. The 
output per architecture is the angular spread for every 
LVLH point at every instant of time, which serves as inputs 
into a science evaluation model to determine how much the 
spread reduces BDRF uncertainties. The constellation 
software is capable of generating angular, spatial and 
temporal metrics, because the analysis is in the global, 
ECEF frame. As with formations, constellation architectures 
can be selected either by comparing the intermediate metrics 
like angular spread or revisit rate or by passing it into 
OSSEs to compare science-based uncertainties. 

There is quantitative dependence between angular and 
temporal sampling (metrics in this paper) and spatial, 
spectral and radiometric resolution (introduced and not 
discussed). Previous studies have discussed this dependence 
for specific instruments in the context of muli-angular 
remote sensing[6]. Future work includes integrating those 
tools with the ones mentioned here, for more exhaustive 
DSM tradespace software. The future tool will also allow 
users to scale up formation clusters to more than one, where 
clusters will provide angular spread and the clustellation 
will provide frequent observations of that spread. 

The described tool has been applied to design and validate 
DSM architectures for BRDF[13], albedo[12] and ERB[11] 
estimation and can be extended to more studies in the future. 


References 

[1] E. M. C. Kong and D. W. Miller, “Optimal spacecraft 
reorientation for earth orbiting clusters: applications to 
Techsat 21,” Acta Astronaut, vol. 53, no. 11, pp. 863- 
877, Dec. 2003. 

[2] S. Nag, C. K. Gatebe, and O. L. De Week, “Relative 
Trajectories for Multi- Angular Earth Observation using 
Science Performance Optimization,” in IEEE Xplore, 
Aerospace Conference 2014, Big Sky, Montana, 2014. 

[3] M. G. O’Neill, H. Yue, S. Nag, P. Grogan, and O. de 
Week, “Comparing and Optimizing the DARPA 
System F6 Program Value-Centric Design 
Methodologies,” in Proceedings of the AIAA Space 
Conference, Anaheim, California, 2010. 

[4] Aleksandr A. Kerzhner, Michel D. Ingham, 
Mohammed O. Khan, Jaime Ramirez, Javier De Luis, 
Jeremy Hollman, Steven Arestie, and David Sternberg, 
“Architecting Cellularized Space Systems using Model- 
Based Design Exploration,” in AIAA SPACE 2013 
Conference and Exposition, 0 vols., American Institute 
of Aeronautics and Astronautics, 2013. 

[5] K. T. Alfriend, S. R. Vadali, and H. Schaub, 

“Formation flying satellites: Control by an 

astrodynamicist,” Celest. Mech. Dyn. Astron., vol. 81, 
no. 1-2, pp. 57-62, 2001. 

[6] S. Nag, K. Cahoy, O. de Week, C. Gatebe, B. Pasquale, 

G. Georgiev, T. Hewagama, and S. Aslam, “Evaluation 
of Hyperspectral Snapshot Imagers onboard 

Nanosatellite Clusters for Multi-Angular Remote 
Sensing,” in AIAA Space Conference, San Diego, 2013. 

[7] S. Nag, “Design of Nano-satellite Cluster Formations 
for Bi-Directional Reflectance Distribution Function 
(BRDF) Estimations,” AIAAUSU Conf Small Satell, 
Aug. 2013. 

[8] D. J. Diner, J. C. Beckert, T. H. Reilly, C. J. Bruegge, J. 
E. Conel, R. A. Kahn, J. V. Martonchik, T. P. 
Ackerman, R. Davies, S. A. W. Gerstl, and others, 
“Multi-angle Imaging SpectroRadiometer (MISR) 
instrument description and experiment overview,” 
Geosci. Remote Sens. IEEE Trans. On, vol. 36, no. 4, 
pp. 1072-1087, 1998. 

[9] X. Xiong, R. Wolfe, W. Barnes, B. Guenther, E. 
Vermote, N. Saleous, and V. Salomonson, “Terra and 
Aqua MODIS Design, Radiometry, and Geometry in 
Support of Land Remote Sensing,” Land Remote Sens. 
Glob. Environ. Change, pp. 133-164, 2011. 

[10] S. Nag, K. Cahoy, and O. L. De Week, “Subsystem 
Support Feasibility for Distributed Space Missions in 
Earth Observation,” in IEEE Xplore, Aerospace 
Conference 2015, Big Sky, Montana, USA, 2015. 

[11] S. Nag, “Satellite Constellation Mission Design using 
Model-Based Systems Engineering and Observing 
System Simulation Experiments,” in Proceedings of the 
Small Satellite Conference, Logan, Utah, 2014. 

[12] S. Nag, C. Gatebe, D. W. Miller, and O. L. De Week, 
“Effect of Satellite Formation Architectures and 
Imaging Modes on Albedo Estimation of major 
Biomes,” Acta Astronaut, 2015. 


[13] S. Nag, C. K. Gatebe, and O. L. De Week, “Observing 
System Simulations for Small Satellite Clusters 
estimating Bi-Directional Reflectance,” J. Appl. 
Geophys., 2015. 

[14] O. L. De Week, U. Scialom, and A. Siddiqi, “Optimal 
reconfiguration of satellite constellations with the 
auction algorithm,” Acta Astronaut., vol. 62, no. 2, pp. 
112-130, 2008. 

[15] S. Nag, J. A. Hoffman, and O. L. de Week, 
“Collaborative and Educational Crowdsourcing of 
Spaceflight Software using SPHERES Zero Robotics,” 
Int. J. Space Technol. Manag. Innov. IJSTMI, vol. 2, 
no. 2, pp. 1-23, 2012. 

[16] S. Nag, “Collaborative competition for crowdsourcing 
spaceflight software and STEM education using 
SPHERES Zero Robotics,” Massachusetts Institute of 
Technology, 2012. 

[17] S. Nag, I. Heffan, A. Saenz-Otero, M. Lydon, 
“SPHERES Zero Robotics software development: 
Lessons on crowdsourcing and collaborative 
competition,” in Aerospace Conference, 2012 IEEE. 

[18] J. Esper, S. Neeck, W. Wiscombe, M. Ryschkewitsch, 
and J. Andary, “Leonardo-BRDF : A New Generation 
Satellite Constellation,” 2000. 

[19] F. E. Nicodemus, Geometrical considerations and 
nomenclature for reflectance, vol. 160. US Department 
of Commerce, National Bureau of Standards 
Washington, D. C, 1977. 

[20] A. Lyapustin, C. K. Gatebe, R. Kahn, et al, “Analysis 
of snow bidirectional reflectance from ARCTAS 
Spring-2008 Campaign,” Atmos Chem Phys, vol. 10, 
no. 9, pp. 4359^375,2010. 

[21] S. Liang, Advances in land remote sensing: System, 
modelling, inversion and application. Springer, 2008. 

[22] B. A. Wielicki and E. F. Harrison, “Mission to planet 
Earth: Role of clouds and radiation in climate,” Bull. 
Am. Meteorol. Soc., vol. 76, no. 11, 1995. 

[23] S. C. Spangelo, D. Kaslow, C. Delp, B. Cole, L. 
Anderson, E. Fosse, B. S. Gilbert, L. Hartman, T. 
Kahn, and J. Cutler, “Applying model based systems 
engineering (MBSE) to a standard CubeSat,” in 
Aerospace Conference, 2012 IEEE, 2012, pp. 1-20. 

[24] P. T. Grogan, “A flexible, modular approach to 
integrated space exploration campaign logistics 
modeling, simulation, and analysis,” Massachusetts 
Institute of Technology, 2010. 

[25] H. Heidt, J. Puig-Suari, A. Moore, S. Nakasuka, and R. 
Twiggs, “CubeSat: A new generation of picosatellite 
for education and industry low-cost space 
experimentation,” 2000. 

[26] O. L. De Week, R. D. Neufville, and M. Chaize, 
“Staged deployment of communications satellite 
constellations in low earth orbit,” J. Aerosp. Comput. 
Inf. Commun., vol. 1, no. 3, pp. 119-136, 2004. 

[27] G. W. Hill, “Researches in the lunar theory,” Am. J. 
Math., vol. 1, no. 1, pp. 5-26, 1878. 

[28] W. H. Clohessy, R.S. Wiltshire, “Terminal Guidance 
System for Satellite Rendezvous,” J. Aerosp. Sci., vol. 
27, pp. 653-658, 1960. 


19 


[29] J. R. Wertz, Orbit & Constellation Design & 
Management, second printing ed. El Segundo. 
California: Microcosm Press, 2009. 

[30] S. A. Schweighart, R. J. Sedwick, “High-Fidelity 
Linearized J Model for Satellite Formation Flight,” J. 
Guid. Control Dyn., vol. 25, no. 6, 2002. 

[31] M. Sabatini and G. B. Palmerini, “Linearized 
formation-flying dynamics in a perturbed orbital 
environment,” in Aerospace Conference, 2008 IEEE, 
2008, pp. 1-13. 

[32] “Systems Tool Kit,” Graphics Inc. [Online]. 

[33] H. Schaub, “Spacecraft relative orbit geometry 
description through orbit element differences,” in 14th 
US National Congress of Theoretical and Applied 
Mechanics, Blacksburg, VA, 2002. 

[34] S. W. Pack, “Reconfigurable satellite constellations for 
geo-spatially adaptive Earth observation missions,” 
Massachusetts Institute of Technology, 2012. 

[35] R.Legge, “Optimization and Valuation of 
Reconfigurable Satellite Constellations Under 
Uncertainty.” MIT Internal Thesis Proposal, 2013. 

[36] T. J. Lang and W. S. Adams, “A comparison of satellite 
constellations for continuous global coverage,” in 
Mission Design & Implementation of Satellite 
Constellations, Springer, 1998, pp. 51-62. 

[37] T. J. Lang, “Walker constellations to minimize revisit 
time in low earth orbit.,” Adv. Astronaut. Sci., vol. 1 14, 
p. 16, 2003. 

[38] 65 Authors from the Astronautics Community, Space 
Mission Engineering: The New SMAD, First. 
Microcosm Press, 2011. 

[39] D. Mortari, M. P. Wilkins, and C. Bruccolerr, “The 
flower constellations.,” Adv. Astronaut. Sci., vol. 115, 
pp. 269-290, 2003. 

[40] M. P. Wilkins, C. Bruccoleri, and D. Mortari, 
“Constellation Design Using Flower Constellations,” 
Pap. AAS, pp. 04-208, 2004. 

[41] M. P. Wilkins and D. Mortari, “Constellation Design 
via Projection of an Arbitrary Shape onto a Flower 
Constellation Surface,” in Paper of the 2004 AIAA/AAS 
Astro dynamics Specialist Conference, Providence, 
Rhode Island, 2004. 

[42] M. E. Avendano and D. Mortari, “New Insights on 
Flower Constellations Theory,” Aerosp. Electron. Syst. 
IEEE Trans. On, vol. 48, no. 2, pp. 1018-1030, 2012. 

[43] S. Nag, J. LeMoigne, and O. L. De Week, “Cost and 
Risk Analysis of Small Satellite Constellations for 
Earth Observation,” in IEEE Xplore, Aerospace 
Conference 2014, Big Sky, Montana, USA, 2014. 

[44] M. KING, M. STRANGE, et al, “Multiwavelength 
scanning radiometer for airborne measurements of 
scattered radiation within clouds,” J. Atmospheric 
Ocean. Technol, vol. 3, pp. 513-522, 1986. 

[45] C. K. Gatebe, “Airborne spectral measurements of 
surface-atmosphere anisotropy for several surfaces and 
ecosystems over southern Africa,” J. Geophys. Res., 
vol. 108, no. D13,2003. 


Biographies 


Sreeja Nag is a PhD candidate in Space 
Systems Engineering at the 
Massachusetts Institute of Technology 
and a part-time research engineer at 
NASA Ames Research Center. She has a 
dual SM in Aeronautics & Astronautics 


Jacqueline Le Moigne is the Assistant 
Chief for Technology in the Software 
Engineering Division at NASA Goddard 
where she is currently leading a study 
on Distributed Spacecraft Missions. She 
has performed extensive work in 
developing new technologies for remote sensing data 
analysis, e.g., image registration, high-performance and 
onboard processing. She has published over 100 
publications and recently co-edited a book on 'Jmage 
Registration for Remote Sensing. ” Email: 
Jacqueline.LeMoigne&nasa.gov 



David W. Miller is a Hunsacker 
Professor of Aeronautics and 
Astronautics at the 

Massachusetts Institute of Technology. 
He is currently serving as the Chief 
Technologist at NASA Headquarters. As 
the Director of the Space Systems Laboratory at MIT, he 
was the PI of the SPHERES program aboard the 
International Space Station and, more recently, served as 
the vice chair of the Air Force Scientific Advisory Board. 
Email: miller d(d)mit. edu 



Olivier de Week is a Professor of 
Aeronautics and Astronautics and 
Engineering Systems at the 
Massachusetts Institute of Technology. 
He is also the Executive Director of 
MIT Production in the Innovation 



Engineering along with Technology & 
Policy at MIT. She has worked at NASA Goddard Space 
Flight Center, has summer research experience with NASA 
JPL in 2008, the European Space Agency (ESTEC) in 2010 
and led the SPHERES Zero Robotics Program in 2011. 
Email: sreeia_n(q)mit. edu 



Economy (PIE) Study and the Co- 
Director, Center for Complex Engineering 
Systems at KACST and MIT. Email: deweck(d),mit. edu 


20 





